# PACKAGES
## Tidy code
library(dplyr)
library(magrittr)
library(tibble)
library(purrr)
## Plotting
library(ggplot2)
library(patchwork)
library(reshape2)
library(ggforce)
library(ggpubr)
library(plotly)
library(graphics)
library(circlize)
## Stats
library(MASS)
library(car)
library(FactoMineR)
library(factoextra)
library(ca)
library(vegan)
library(cluster)
library(heplots)
library(rstatix)
library(MVN)
# custom functions
profile_dist_matrix <- function(data){
mat <- as.matrix((data))
n <- ncol(mat)
profiles <- prop.table(x = mat, margin = 2)
average.profile <- rowSums(mat)/sum(mat)
dist.mat<- matrix(NA, n, n)
diag(dist.mat) <- 0
for (i in 1:(n - 1)) {
for (j in (i + 1):n) {
d2 <- sum(((profiles[, i] - profiles[, j])^2) / average.profile)
dist.mat[i, j] <- dist.mat[j, i] <- d2
}
}
colnames(dist.mat) <- rownames(dist.mat) <- colnames(mat)
dist.mat
}
El código y los archivos orginales de esta PEC, así como los informes en formato HTML pueden consultarse en el siguiente repositorio de GitHub. El repositorio es privado, así que se tendrá que requerir el acceso para poder ver el código.
En el trabajo de Hunt et al. se estudió la capacidad reproductiva de cinco especies de aves marinas en dos colonias en el sureste del mar de Bering. Además, el apéndice de este estudio resume las colonias y los tamaños de las poblaciones de otros trabajos. El archivo seabirds.csv recoge los datos (número de pájaros) de 23 especies en 9 colonias en el área del norte polar y subpolar. El principal interés de este ejercicio es representar las colonias de diversas formas y estudiar posibles conglomerados.
seabirds <- read.csv("seabirds.csv", stringsAsFactors = F)
seabirds2 <- seabirds %>% column_to_rownames("Specie")
margin_col <- rowSums(seabirds2)
margin_row <- colSums(seabirds2)
n <- sum(seabirds2)
Calcular las frecuencias relativas, las frecuencias relativas marginales y la matriz de perfiles. El resultado debería ser la tabla 12.6 del libro de Krebs y que reproducimos al final de este documento.
Para calcular las frecuencias relativas, solo tenemos que dividir cada valor (frecuencias absolutas) por el total de observaciones (4765889) o usando la función prop.table(x=seabirds) y las frecuencias relativas marginales se consiguen con la suma de las frecuencias relativas de cada columna o de cada fila. Sin embargo, con esto no nos sale igual que la tabla 12.6 del libro de Krebs.
seabirds_freq_rel2 <- prop.table(x = seabirds2 %>% as.matrix())
addmargins(seabirds_freq_rel2) %>% round(4) %>%
knitr::kable(caption = "Table 1.1. Relative and marginal frequencies from data in 'seabirds.csv'. Relative frequencies are calculated using the total number of observations. Marginal frequencies are observed in the last row and the last column with the column/row name 'sum'. This table is not coincident with table 12.6 from Krebs book.", align = "c")
| CH | PLI | CI | NS | CL | CT | SI | SPI | SGI | Sum | |
|---|---|---|---|---|---|---|---|---|---|---|
| Northern.fulmar | 0.0000 | 0.0260 | 0.0000 | 0.0000 | 0.0000 | 0.0000 | 0.0000 | 0.0001 | 0.0147 | 0.0409 |
| Glaucous-winged.gull | 0.0000 | 0.0001 | 0.0000 | 0.0001 | 0.0000 | 0.0001 | 0.0000 | 0.0000 | 0.0000 | 0.0003 |
| Black-legged.kittiwake | 0.0084 | 0.0122 | 0.0126 | 0.0016 | 0.0052 | 0.0056 | 0.0010 | 0.0065 | 0.0151 | 0.0682 |
| Red-legged.kittiwake | 0.0000 | 0.0000 | 0.0000 | 0.0000 | 0.0000 | 0.0000 | 0.0000 | 0.0005 | 0.0462 | 0.0466 |
| Thick-billed.murre | 0.0588 | 0.0361 | 0.0671 | 0.0001 | 0.0063 | 0.0490 | 0.0000 | 0.0231 | 0.3147 | 0.5552 |
| Common.murre | 0.0000 | 0.0000 | 0.0000 | 0.0088 | 0.0147 | 0.0327 | 0.0011 | 0.0082 | 0.0399 | 0.1053 |
| Black.guillemot | 0.0000 | 0.0017 | 0.0000 | 0.0000 | 0.0000 | 0.0000 | 0.0000 | 0.0000 | 0.0000 | 0.0018 |
| Pigeon.guillemot | 0.0000 | 0.0000 | 0.0000 | 0.0000 | 0.0000 | 0.0000 | 0.0000 | 0.0000 | 0.0000 | 0.0000 |
| Horned.puffin | 0.0000 | 0.0000 | 0.0000 | 0.0007 | 0.0003 | 0.0003 | 0.0000 | 0.0009 | 0.0059 | 0.0081 |
| Tufted.puffin | 0.0000 | 0.0000 | 0.0000 | 0.0000 | 0.0000 | 0.0000 | 0.0000 | 0.0002 | 0.0013 | 0.0015 |
| Atlantic.puffin | 0.0000 | 0.0000 | 0.0000 | 0.0000 | 0.0000 | 0.0000 | 0.0034 | 0.0000 | 0.0000 | 0.0034 |
| Pelagic.cormorant | 0.0000 | 0.0000 | 0.0000 | 0.0001 | 0.0000 | 0.0000 | 0.0000 | 0.0000 | 0.0000 | 0.0001 |
| Red-faced.cormorant | 0.0000 | 0.0000 | 0.0000 | 0.0000 | 0.0000 | 0.0000 | 0.0000 | 0.0005 | 0.0010 | 0.0016 |
| Shag | 0.0000 | 0.0000 | 0.0000 | 0.0000 | 0.0000 | 0.0000 | 0.0000 | 0.0000 | 0.0000 | 0.0000 |
| Parakeet.auklet | 0.0000 | 0.0000 | 0.0000 | 0.0000 | 0.0000 | 0.0000 | 0.0000 | 0.0071 | 0.0315 | 0.0386 |
| Crested.auklet | 0.0000 | 0.0000 | 0.0000 | 0.0000 | 0.0000 | 0.0000 | 0.0000 | 0.0013 | 0.0059 | 0.0071 |
| Least.auklet | 0.0000 | 0.0000 | 0.0000 | 0.0000 | 0.0000 | 0.0000 | 0.0000 | 0.0048 | 0.0525 | 0.0573 |
| Razorbill | 0.0000 | 0.0000 | 0.0000 | 0.0000 | 0.0000 | 0.0000 | 0.0009 | 0.0000 | 0.0000 | 0.0009 |
| Manx.shearwater | 0.0000 | 0.0000 | 0.0000 | 0.0000 | 0.0000 | 0.0000 | 0.0546 | 0.0000 | 0.0000 | 0.0546 |
| Storm.petrel | 0.0000 | 0.0000 | 0.0000 | 0.0000 | 0.0000 | 0.0000 | 0.0027 | 0.0000 | 0.0000 | 0.0027 |
| Herring.gull | 0.0000 | 0.0000 | 0.0000 | 0.0000 | 0.0000 | 0.0000 | 0.0016 | 0.0000 | 0.0000 | 0.0016 |
| Great.black-backed.gull | 0.0000 | 0.0000 | 0.0000 | 0.0000 | 0.0000 | 0.0000 | 0.0001 | 0.0000 | 0.0000 | 0.0001 |
| Lesser.black.backed.gull | 0.0000 | 0.0000 | 0.0000 | 0.0000 | 0.0000 | 0.0000 | 0.0042 | 0.0000 | 0.0000 | 0.0042 |
| Sum | 0.0672 | 0.0760 | 0.0798 | 0.0113 | 0.0266 | 0.0876 | 0.0696 | 0.0533 | 0.5285 | 1.0000 |
Para que coincida con la tabla 12.6, la tabla de frecuencias relativas debe hacerse respecto al total de observaciones en cada lugar, por lo que el valor de la frecuencia absoluta dividido por la suma de todas las observaciones en la columna correspondiente, algo que también pude conseguirse con prop.table(x = seabirds, margin = 2). Lo que es la matriz de perfiles de las columnas o la frecuencia relativa condicionada a las columnas (frec. absoluta de una celda dividida por la suma de las frec. absolutas de la columna correspondiente).
seabirds_perfiles_c <- prop.table(x = seabirds2 %>% as.matrix(), margin = 2)
addmargins(seabirds_perfiles_c) %>% round(4) %>%
knitr::kable(caption = "Table 1.2. Relative and marginal frequencies conditioned to columns from data in 'seabirds.csv'. Relative frequencies are calculated taking all the observations in each place (column), not the total number of observations. Marginal frequencies are observed in the last row and the last column with the column/row name 'sum'.", align = "c")
| CH | PLI | CI | NS | CL | CT | SI | SPI | SGI | Sum | |
|---|---|---|---|---|---|---|---|---|---|---|
| Northern.fulmar | 0.0000 | 0.3422 | 0.0000 | 0.0000 | 0.0000 | 0.0000 | 0.0007 | 0.0028 | 0.0278 | 0.3734 |
| Glaucous-winged.gull | 0.0005 | 0.0011 | 0.0004 | 0.0051 | 0.0004 | 0.0007 | 0.0000 | 0.0000 | 0.0000 | 0.0082 |
| Black-legged.kittiwake | 0.1249 | 0.1600 | 0.1577 | 0.1402 | 0.1972 | 0.0634 | 0.0151 | 0.1221 | 0.0286 | 1.0094 |
| Red-legged.kittiwake | 0.0000 | 0.0000 | 0.0000 | 0.0000 | 0.0000 | 0.0000 | 0.0000 | 0.0087 | 0.0873 | 0.0960 |
| Thick-billed.murre | 0.8740 | 0.4746 | 0.8413 | 0.0074 | 0.2367 | 0.5593 | 0.0000 | 0.4334 | 0.5955 | 4.0222 |
| Common.murre | 0.0000 | 0.0000 | 0.0000 | 0.7765 | 0.5522 | 0.3728 | 0.0160 | 0.1537 | 0.0754 | 1.9466 |
| Black.guillemot | 0.0006 | 0.0221 | 0.0005 | 0.0000 | 0.0013 | 0.0000 | 0.0000 | 0.0000 | 0.0000 | 0.0246 |
| Pigeon.guillemot | 0.0000 | 0.0000 | 0.0000 | 0.0000 | 0.0000 | 0.0000 | 0.0000 | 0.0000 | 0.0000 | 0.0000 |
| Horned.puffin | 0.0000 | 0.0000 | 0.0000 | 0.0592 | 0.0114 | 0.0036 | 0.0000 | 0.0173 | 0.0111 | 0.1026 |
| Tufted.puffin | 0.0000 | 0.0000 | 0.0000 | 0.0008 | 0.0002 | 0.0000 | 0.0000 | 0.0039 | 0.0024 | 0.0072 |
| Atlantic.puffin | 0.0000 | 0.0000 | 0.0000 | 0.0000 | 0.0000 | 0.0000 | 0.0482 | 0.0000 | 0.0000 | 0.0482 |
| Pelagic.cormorant | 0.0000 | 0.0000 | 0.0000 | 0.0096 | 0.0006 | 0.0001 | 0.0001 | 0.0000 | 0.0000 | 0.0104 |
| Red-faced.cormorant | 0.0000 | 0.0000 | 0.0000 | 0.0000 | 0.0000 | 0.0000 | 0.0000 | 0.0099 | 0.0020 | 0.0118 |
| Shag | 0.0000 | 0.0000 | 0.0000 | 0.0000 | 0.0000 | 0.0000 | 0.0001 | 0.0000 | 0.0000 | 0.0001 |
| Parakeet.auklet | 0.0000 | 0.0000 | 0.0000 | 0.0012 | 0.0000 | 0.0000 | 0.0000 | 0.1340 | 0.0595 | 0.1947 |
| Crested.auklet | 0.0000 | 0.0000 | 0.0000 | 0.0000 | 0.0000 | 0.0000 | 0.0000 | 0.0236 | 0.0111 | 0.0348 |
| Least.auklet | 0.0000 | 0.0000 | 0.0000 | 0.0000 | 0.0000 | 0.0000 | 0.0000 | 0.0906 | 0.0992 | 0.1899 |
| Razorbill | 0.0000 | 0.0000 | 0.0000 | 0.0000 | 0.0000 | 0.0000 | 0.0130 | 0.0000 | 0.0000 | 0.0130 |
| Manx.shearwater | 0.0000 | 0.0000 | 0.0000 | 0.0000 | 0.0000 | 0.0000 | 0.7838 | 0.0000 | 0.0000 | 0.7838 |
| Storm.petrel | 0.0000 | 0.0000 | 0.0000 | 0.0000 | 0.0000 | 0.0000 | 0.0389 | 0.0000 | 0.0000 | 0.0389 |
| Herring.gull | 0.0000 | 0.0000 | 0.0000 | 0.0000 | 0.0000 | 0.0000 | 0.0229 | 0.0000 | 0.0000 | 0.0229 |
| Great.black-backed.gull | 0.0000 | 0.0000 | 0.0000 | 0.0000 | 0.0000 | 0.0000 | 0.0009 | 0.0000 | 0.0000 | 0.0009 |
| Lesser.black.backed.gull | 0.0000 | 0.0000 | 0.0000 | 0.0000 | 0.0000 | 0.0000 | 0.0603 | 0.0000 | 0.0000 | 0.0603 |
| Sum | 1.0000 | 1.0000 | 1.0000 | 1.0000 | 1.0000 | 1.0000 | 1.0000 | 1.0000 | 1.0000 | 9.0000 |
Calcular la matriz de distancias ji-cuadrado entre los perfiles de las columnas y su inercia total.
Los perfiles de las columnas estan calculados en el ejercicio enterior, así que usaremos esos datos para calcular las distancias ji-cuadrado y su inercia total.
Distancia ji-cuadrado: con la distancia ji-cuadrado entre los perfiles, podemos observar qué columnas se parecen más unas a otras: a una distancia más pequeña, más similitud entre columnas (en este caso, colonias de varias aves).
La distancia ji-cuadrado entre columnas se puede definir como
\[ d_{ij}^{(cols)} = \sum_{k=1}^{r} \frac{1}{p_k.}(p^c_{ki} - p^c_{kj})^2 \]
donde:
En nuestro caso, las proporciones las tenemos en la tabla calculada en el ejercicio anterior (ejercicio 1a) y que coinciden con la tabla 12.6 del libro de Krebs.
Para calcular dicha matriz, usaremos la función profile_dist_matrix(), definida a l’inicio de este informe, sacada de STHDA y modificada por mi. profile_dist_matrix() acepta como input la matriz de perfiles (que devemos de transformar con t() para que coja las columnas) y el perfil de columnas medio (\(p_k.\) en la función de la distancia entre columnas)
dist_cols_matrix <- profile_dist_matrix((seabirds2))
dist_cols_matrix %>% knitr::kable(caption = "Table 1.3. Chi-squared distance matrix between column profiles of the 'seabirds.csv' data.")
| CH | PLI | CI | NS | CL | CT | SI | SPI | SGI | |
|---|---|---|---|---|---|---|---|---|---|
| CH | 0.0000000 | 3.425238 | 0.0178094 | 8.251115 | 3.7234212 | 1.5564260 | 15.46510 | 1.3725624 | 0.8159060 |
| PLI | 3.4252376 | 0.000000 | 3.3647355 | 10.407262 | 6.1401304 | 4.6055964 | 17.75168 | 4.1374664 | 3.4917128 |
| CI | 0.0178094 | 3.364736 | 0.0000000 | 8.154352 | 3.5963874 | 1.5964295 | 15.48546 | 1.3407506 | 0.8934179 |
| NS | 8.2511147 | 10.407262 | 8.1543518 | 0.000000 | 1.5595341 | 3.2777270 | 20.80417 | 5.7319920 | 6.9746456 |
| CL | 3.7234212 | 6.140130 | 3.5963874 | 1.559534 | 0.0000000 | 0.7664893 | 17.24585 | 2.4287149 | 3.2811626 |
| CT | 1.5564260 | 4.605596 | 1.5964295 | 3.277727 | 0.7664893 | 0.0000000 | 15.71938 | 1.3213094 | 1.3394868 |
| SI | 15.4650970 | 17.751675 | 15.4854630 | 20.804173 | 17.2458462 | 15.7193804 | 0.00000 | 15.3926081 | 15.0680100 |
| SPI | 1.3725624 | 4.137466 | 1.3407506 | 5.731992 | 2.4287149 | 1.3213094 | 15.39261 | 0.0000000 | 0.5942528 |
| SGI | 0.8159060 | 3.491713 | 0.8934179 | 6.974646 | 3.2811626 | 1.3394868 | 15.06801 | 0.5942528 | 0.0000000 |
Como és de esperar, en la diagonal observamos distancias de 0, ya que comparamos cada columna consigo misma.
Por otro lado, podemos ver como hay columnas con una distancia muy pequeña como CH-CI (0.0178) o CL-CT (0.7665), mientras que otras tienen una distancia mayor y, por lo tanto, observamos como són distintas (CH-SI = 15.4651, NS-SI = 20.8042).
Inercia total: la inercia total es la media ponderada de los cuadrados de las distancias \(\chi^2\) entre los perfiles fila/columna y el perfil media [^(2)^](https://www.fbbva.es/wp-content/uploads/2017/05/dat/greenacre_cap04.pdf). También puede encontrarse dividiendo el estadístico \(\chi^2\) por el total de observaciones. Cabe destacar que la inercia total es la misma tanto si se hace con los perfiles de fila como los de columna, ya que estamos actuando sobre la totalidad de los datos.
as.numeric(chisq.test(seabirds2)$stat)/n
## [1] 1.565692
Con la matriz de distancias ji-cuadrado entre los perfiles realizar un escalado multidimensional. Dibujar las coordenadas principales para las columnas.
Para el escalado multidimensional usaremos la función cmdscale() con 8 coordenadas principales (el máximo en este caso) y indicando que queremos que nos devuelva los valores propios (eigenvalues).
c <- cmdscale(d = dist_cols_matrix, k=8, eig = T)
c
## $points
## [,1] [,2] [,3] [,4] [,5]
## CH 0.5043776 -1.75638463 -0.72114533 0.02370077 0.120329377
## PLI 1.9436165 -4.95136344 0.48741897 0.10467367 -0.010646143
## CI 0.5461064 -1.70615994 -0.70169181 -0.09456648 -0.110423324
## NS 6.0813689 4.66506492 0.15933028 0.03698951 0.001837342
## CL 2.8046648 1.78390484 -0.31619954 -0.11640599 -0.016327900
## CT 1.4177905 1.07090468 0.05751178 0.48401878 -0.004345404
## SI -14.5809325 1.59022765 0.06534089 0.02022302 -0.002890291
## SPI 0.8713107 0.05057765 0.41238536 -0.57607534 0.021088565
## SGI 0.4116971 -0.74677172 0.55704940 0.11744207 0.001377778
##
## $eig
## [1] 2.647219e+02 5.969290e+01 1.863310e+00 6.157183e-01 2.752965e-02
## [6] -3.197442e-14 -6.373612e-01 -1.598085e+00 -1.998594e+01
##
## $x
## NULL
##
## $ac
## [1] 0
##
## $GOF
## [1] 0.9363544 1.0000000
Una vez realizado el escalado multidimensional, debemos ver cuantas dimensiones (o coordenadas principales) explican la mayor parte de los datos originales.
Para ello usaremos los dos métodos siguientes, que se basan en la suma acumulada de los valores propios dividida por la suma total de estos valores propios. Como en nuestro caso hay algunos valores negativos, es necesario usar el valor absoluto (método 1) o elevar al cuadrado (método 2).
\[ P_{m}^{(1)} = \frac{\sum^{m}_{i=1} |\lambda_{i}| }{\sum^{n}_{i=1} |\lambda_{i}| } \]
cumsum(abs(c$eig)) / sum(abs(c$eig))
## [1] 0.7582053 0.9291753 0.9345121 0.9362756 0.9363544 0.9363544 0.9381799
## [8] 0.9427571 1.0000000
\[ P_{m}^{(1)} = \frac{\sum^{m}_{i=1} \lambda_{i}^2 }{\sum^{n}_{i=1} \lambda_{i}^2 } \]
cumsum(c$eig^2) / sum(c$eig^2)
## [1] 0.9463924 0.9945136 0.9945605 0.9945656 0.9945657 0.9945657 0.9945711
## [8] 0.9946056 1.0000000
Ambos métodos nos sugieren utilizar dos dimensiones para la representación de los datos, pues después de la segunda componente principal, las otras no parecen tener un gran papel en la explicación de los datos
c$points %>% as.data.frame() %>% mutate(names = rownames(c$points)) %>%
ggplot(aes(V1, V2)) +
geom_point() +
geom_text(aes(label = names),
position = position_dodge(width = 1),
vjust = -1, hjust = 1, size = 4) +
ggtitle("MDS of Chi-squared distances",
"of column profiles") +
xlab("1st principal coordinate") + ylab("2nd principal coordinate") +
theme_bw(base_size = 15)
Figure 1.1. MDS of chi-squared distances of column profiles. 2D representation (dimensións 1 and 2).
Como se puede ver, la mayoría de columnas se encuentran en un punto similar dentro de la primera coordinada principal (salvo la columna SI), que es la que explica la mayor parte de la variación, mientras que en la segunda coordinada varian mucho.
Vemos como con dos dimensiones se explica mucho la distancia/dissimilaridad entre columnas, ya que cuanto más alejadas están unas de otras en el gráfico de coordenadas principales, más distancia hay entre sus perfiles. Por ejemplo, SI vemos como tiene un perfil que tiene una distáncia grande con todas las otras, PLI y NS también se encuentran lejos entre sí, mientras que CI y CH están muy cerca.
A continuación, aunque no es extrictamente necesario, vamos ha hacer un gráfico 3D con las primeras tres coordinadas principales.
plot3d <- cmdscale(dist_cols_matrix, k=8) %>% as.data.frame() %>% rownames_to_column("Name")
plotly::plot_ly(x = plot3d$V1, y = plot3d$V2, z = plot3d$V3, name = plot3d$Name,
type = "scatter3d")
Figura 1.2. MDS of chi-squared distances of column profiles. 3D representation (dimensións 1, 2 and 3).
Podemos ver como en el eje Z (tercera coordinada principal), las columnas no parecen tener mucha variabilidad, dado que el rango de Z es de 1, mientras que los otros dos son de 20 y 8 para X e Y respectivamente.
Aún así, vemos como esta nueva coordenada principal explica un poco más la distancia entre los perfiles de las columnas. Por ejemplo, las columnas CI y CH, que con dos dimensiones se veían muy cerca de SGI, ahora se encuentran en los extremos opuestos del eje Z.
Sin embargo, esta representación puede ser confusa, ya que debido a las escalas del gráficovemos como, por ejemplo, CH se encuentra muy distanciada de SGI y parece estar más cerca de CL, cuando las distancias ji-cuadrado entre sus perfiles indican lo contrario.
Realizar un análisis de correspondencias y calcular las inercias principales (en %) y la inercia total con los valores propios. Dibujar una representación simétrica del CA. A pesar de la confusión de nombres, ¿cuales son las especies que caracterizan a la colonia SI (Skomer Island, Irish Sea)?
El análisis de correspondencias vamos a hacerlo usando la función CA() del paquete FactoMineR. Esta función nos devuelve los valores propios (eigenvalues) y los porcentajes de variancia (inercia) de cada dimensión ademas de las coordenadas de las filas y las columnas en cada dimensión y sus inercias, contribuciones, etc. Asimismo, nos permite hacer directamente una representación de los datos (indicando graph=T).
Como el porcentaje de inercia ya nos lo da la función, no hace falta calcularlo, aunque sencillamente sería calcular la suma de los valores propios (inercia total) y efectuar el porcentaje de cada valor propio respecto a ese valor.
\[ Inercia \space total = {\sum_{i=1}^{n}\lambda_i}*100 \] Siendo \(\lambda_i\) los valores propios y \(n\) el total de dimensiones del anàlisis de correspondencias.
\[ Inercia \space principal \space _{(de \space la \space dimension \space j)} = \frac{\lambda_j}{\sum_{i=1}^n\lambda_i}*100 \]
ca = CA(X = seabirds2, ncp = 5, graph = F)
Los valores propios, así como la inercia principal (en porcentaje) y el porcentaje cumulativo de variancia/inercia se encuentran dentro la el resultado de la funcón CA(), concretamente en la variable $eig.
ca$eig %>% knitr::kable(caption = "Table 1.4. Eigenvalues, percentages of variance/inertia and cumulative percentage of variance/inertia for each dimention.")
| eigenvalue | percentage of variance | cumulative percentage of variance | |
|---|---|---|---|
| dim 1 | 0.9662498 | 61.7139287 | 61.71393 |
| dim 2 | 0.2566386 | 16.3913870 | 78.10532 |
| dim 3 | 0.2059249 | 13.1523299 | 91.25765 |
| dim 4 | 0.0986173 | 6.2986387 | 97.55628 |
| dim 5 | 0.0261770 | 1.6719151 | 99.22820 |
| dim 6 | 0.0083129 | 0.5309395 | 99.75914 |
| dim 7 | 0.0037650 | 0.2404718 | 99.99961 |
| dim 8 | 0.0000061 | 0.0003892 | 100.00000 |
En los datos anteriores podemos ver los distintos porcentajes de inercia de cada dimensión y se ve claramente que la dimensión 1 nos explica la mayor parte de la variación (\(\pm\) 61.7%), mientras que las dimensiones 2 y 3 (\(\pm\) 16.4 y \(\pm\) 13.2, respectivamente) nos explican un porcentaje relativamente mayor que las siguientes ddimensiones, que nos explican muy poco en total.
Para calcular la inercia total tenemos que sumar todos los valores propios, tal y como hemos mencionado anteriormente.
paste("La inercia total es:", sum(ca$eig[,1]))
## [1] "La inercia total es: 1.56569159050875"
Finalmente, para representar los datos podemos usar las coordenadas de las filas y las columnas en las dos primeras dimensiones y hacer un plot() de ellas o simplemente indicando graph=T dentro de la función CA() o usando ggplot2.
ggplot() +
geom_point(ca$row$coord %>% as.data.frame(), mapping = aes(x = `Dim 1`, y = `Dim 2`), color="Blue") +
ggrepel::geom_text_repel(ca$row$coord %>% as.data.frame() %>% rownames_to_column("Name"),
mapping = aes(x = `Dim 1`, y = `Dim 2`, label = Name), color = "Darkblue", seed = 3) +
geom_point(ca$col$coord %>% as.data.frame(), mapping = aes(x = `Dim 1`, y = `Dim 2`), color="Red") +
ggrepel::geom_label_repel(ca$col$coord %>% as.data.frame() %>% rownames_to_column("Name"),
mapping = aes(x = `Dim 1`, y = `Dim 2`, label = Name), color = "Darkred", seed = 3) +
theme_bw(base_size = 14) +
ggtitle("Representation of correspondence analysis",
"Dimension 1 and dimension 2") +
xlab("Dimension 1 (61.7% of variance)") +
ylab("Dimension 2 (16.4% of variance)")
Figure 1.3. CA representation (dimensions 1 and 2) from ‘seabirds.csv’, with rows (species) represented in blue and columns (colonies) represented in red.
Como se puede ver en la figura 1.3, hay muchas especies características de la colonia SI, entre las cuales se encuentran:
Dada la gran cantidad de ceros en la tabla 12.6, en el libro de Krebs se sugiere la utilización de la distancia de Canberra entre las columnas de la tabla 12.6. La distancia de Canberra no tiene una única definición y, además, ha cambiado a lo largo de la historia. Una posible definición entre dos vectores \(\textbf{p} = (p_1, p_2, . . . , p_k)'\) y \(\textbf{q} = (q_1, q_2, . . . , q_k)'\) de la misma longitud es
\[ d_C(\textbf{p},\textbf{q}) = \sum_{i=1}^{k} \frac{\mid p_i - q_i \mid}{|p_i|+|q_i|} \]
Cuando el denominador es cero, el cociente es NaN, y el sumando se elimina.
Comprobar que esta definición no sirve para calcular la matriz de similaridades de la tabla 12.7 del libro de Krebs y que se reproduce al final de este documento.
mat <- as.matrix(seabirds2)
n <- ncol(mat)
profiles <- prop.table(x = mat, margin = 2)
dist.mat<- matrix(NA, n, n)
diag(dist.mat) <- 0
for(i in 1:(n-1)){
for(j in (i+1):n){
dc <- sum(abs(profiles[, i] - profiles[, j]) /
(abs(profiles[, i]) + abs(profiles[, j])), na.rm = T)
dist.mat[i, j] <- dist.mat[j, i] <- dc
}
}
colnames(dist.mat) <- rownames(dist.mat) <- colnames(mat)
dist.mat %>% knitr::kable(caption = "Table 1.5. Canberra distance between column profiles using the 1st definition (mentioned above): d(p,q) = sum(abs(p-q) / (abs(p)+abs(q))), where p and q are 2 column profiles.")
| CH | PLI | CI | NS | CL | CT | SI | SPI | SGI | |
|---|---|---|---|---|---|---|---|---|---|
| CH | 0.0000000 | 2.768917 | 0.3066923 | 7.872672 | 5.248740 | 5.726277 | 14.78458 | 11.34803 | 11.81702 |
| PLI | 2.7689165 | 0.000000 | 2.7128716 | 8.678678 | 6.797518 | 6.724761 | 14.82392 | 11.16369 | 11.65965 |
| CI | 0.3066923 | 2.712872 | 0.0000000 | 7.897405 | 5.108834 | 5.883052 | 14.82557 | 11.44721 | 11.86432 |
| NS | 7.8726718 | 8.678678 | 7.8974045 | 0.000000 | 6.345527 | 8.318004 | 16.74805 | 10.90990 | 11.61912 |
| CL | 5.2487396 | 6.797518 | 5.1088341 | 6.345527 | 0.000000 | 5.609152 | 16.56029 | 11.22128 | 11.82774 |
| CT | 5.7262773 | 6.724761 | 5.8830525 | 8.318004 | 5.609152 | 0.000000 | 15.66589 | 12.51723 | 12.58676 |
| SI | 14.7845812 | 14.823919 | 14.8255729 | 16.748052 | 16.560290 | 15.665893 | 0.00000 | 19.19796 | 18.91136 |
| SPI | 11.3480265 | 11.163693 | 11.4472135 | 10.909901 | 11.221283 | 12.517229 | 19.19796 | 0.00000 | 4.67864 |
| SGI | 11.8170207 | 11.659654 | 11.8643193 | 11.619124 | 11.827739 | 12.586763 | 18.91136 | 4.67864 | 0.00000 |
La tabla 12.7 del libro de Krebs muestra las similaridades entre columnas que se han obtenido calculando \(1 - Distancia\space de \space Canberra\). Como vemos en la tabla 1.5, esta muestra distancias mayores que 1, por lo que restarlas a 1 resultaría en similaridades negativas, cosa que no se muestra en la tabla 12.7 del libro de Krebs.
Una modificación de la distancia anterior es considerar la distancia
\[ d_C(\textbf{p},\textbf{q}) = \frac{1}{k} \sum_{i=1}^{k} \frac{\mid p_i - q_i \mid}{|p_i|+|q_i|} \]
Igual que antes, cuando el denominador es cero, el sumando se elimina.
Comprobar que con esta definición se puede obtener la tabla 12.7.
mat <- as.matrix(seabirds2)
n <- ncol(mat)
k <- nrow(mat)
profiles <- prop.table(x = mat, margin = 2)
dist.mat<- matrix(NA, n, n)
diag(dist.mat) <- 0
for(i in 1:(n-1)){
for(j in (i+1):n){
dc <- sum(abs(profiles[, i] - profiles[, j]) /
(abs(profiles[, i]) + abs(profiles[, j])), na.rm = T) / k
dist.mat[i, j] <- dist.mat[j, i] <- dc
}
}
colnames(dist.mat) <- rownames(dist.mat) <- colnames(mat)
dist.mat %>% knitr::kable(caption = "Table 1.6. Canberra distance between column profiles using the 2nd definition (mentioned above): d(p,q) = (1/k)*(sum(abs(p-q) / (abs(p)+abs(q)))), where p and q are 2 column profiles and k is the number of rows in the profiles table.")
| CH | PLI | CI | NS | CL | CT | SI | SPI | SGI | |
|---|---|---|---|---|---|---|---|---|---|
| CH | 0.0000000 | 0.1203877 | 0.0133344 | 0.3422901 | 0.2282061 | 0.2489686 | 0.6428079 | 0.4933925 | 0.5137835 |
| PLI | 0.1203877 | 0.0000000 | 0.1179509 | 0.3773338 | 0.2955442 | 0.2923809 | 0.6445182 | 0.4853779 | 0.5069415 |
| CI | 0.0133344 | 0.1179509 | 0.0000000 | 0.3433654 | 0.2221232 | 0.2557849 | 0.6445901 | 0.4977049 | 0.5158400 |
| NS | 0.3422901 | 0.3773338 | 0.3433654 | 0.0000000 | 0.2758925 | 0.3616524 | 0.7281762 | 0.4743435 | 0.5051793 |
| CL | 0.2282061 | 0.2955442 | 0.2221232 | 0.2758925 | 0.0000000 | 0.2438762 | 0.7200126 | 0.4878819 | 0.5142495 |
| CT | 0.2489686 | 0.2923809 | 0.2557849 | 0.3616524 | 0.2438762 | 0.0000000 | 0.6811258 | 0.5442274 | 0.5472506 |
| SI | 0.6428079 | 0.6445182 | 0.6445901 | 0.7281762 | 0.7200126 | 0.6811258 | 0.0000000 | 0.8346939 | 0.8222330 |
| SPI | 0.4933925 | 0.4853779 | 0.4977049 | 0.4743435 | 0.4878819 | 0.5442274 | 0.8346939 | 0.0000000 | 0.2034191 |
| SGI | 0.5137835 | 0.5069415 | 0.5158400 | 0.5051793 | 0.5142495 | 0.5472506 | 0.8222330 | 0.2034191 | 0.0000000 |
La tabla 1.6 no muestra ninguna distancia mayor que uno, por lo que es posible que, si restamos estas distancias a 1, el resultado sea la tabla 12.7 del libro de Krebs.
round(1-dist.mat, 2) %>% knitr::kable(caption = "Table 1.7. Table of similarities obtained by resting the Canberra distances in table 1.6 from 1. As it can be seen, this table is the same as the table 12.7 from the Krebs' book.")
| CH | PLI | CI | NS | CL | CT | SI | SPI | SGI | |
|---|---|---|---|---|---|---|---|---|---|
| CH | 1.00 | 0.88 | 0.99 | 0.66 | 0.77 | 0.75 | 0.36 | 0.51 | 0.49 |
| PLI | 0.88 | 1.00 | 0.88 | 0.62 | 0.70 | 0.71 | 0.36 | 0.51 | 0.49 |
| CI | 0.99 | 0.88 | 1.00 | 0.66 | 0.78 | 0.74 | 0.36 | 0.50 | 0.48 |
| NS | 0.66 | 0.62 | 0.66 | 1.00 | 0.72 | 0.64 | 0.27 | 0.53 | 0.49 |
| CL | 0.77 | 0.70 | 0.78 | 0.72 | 1.00 | 0.76 | 0.28 | 0.51 | 0.49 |
| CT | 0.75 | 0.71 | 0.74 | 0.64 | 0.76 | 1.00 | 0.32 | 0.46 | 0.45 |
| SI | 0.36 | 0.36 | 0.36 | 0.27 | 0.28 | 0.32 | 1.00 | 0.17 | 0.18 |
| SPI | 0.51 | 0.51 | 0.50 | 0.53 | 0.51 | 0.46 | 0.17 | 1.00 | 0.80 |
| SGI | 0.49 | 0.49 | 0.48 | 0.49 | 0.49 | 0.45 | 0.18 | 0.80 | 1.00 |
Como se puede ver en la tabla 1.7, las similaridades son las mismas que en la tabla 12.7 del libro de Krebs (aunque la tabla 1.7 es simetrica respecto la diagonal y la del libro de krebs solo muestra la parte superior de esta diagonal). Por lo tanto, podemos afirmar que la definición usada en la tabla 1.6 para calcular las distancias de Canberra es la misma que la que usan los autores de la 12.7 del libro de Krebs.
Finalmente, se puede comprobar que ésta tampoco es la definición que utiliza R para calcular la distancia de Canberra. Tras una ardua investigación, se comprueba que la definición de R es
\[ d_C(\textbf{p},\textbf{q}) = \frac{k}{k-n_z} \sum_{i=1}^{k} \frac{\mid p_i - q_i \mid}{|p_i|+|q_i|} \]
donde nz es el número de denominadores cero.
Comprobar que ésta es la definición de la distancia de Canberra según R.
Según la documentación de la función dist(), el cálculo de la distancia de Canberra usando la función dist(x, method = "canberra") del paquete stats no es la que se menciona arriba, sinó la primera (vease la [documentación dedist()`](https://www.rdocumentation.org/packages/stats/versions/3.6.2/topics/dist)).
Sin embargo, si usamos la definición mencionada en este ejercicio 1.e.3 y calculamos la distancia de Canberra (tabla 1.8) obtenemos las mismas distancias que usando la función dist() de R (tabla 1.9), por lo que podemos afirmar que R usa la definicion mencionada en este ejercicio 1.e.3.
mat <- as.matrix(seabirds2)
n <- ncol(mat)
k <- nrow(mat)
#nz <- sum(mat==0)
profiles <- prop.table(x = mat, margin = 2)
dist.mat<- matrix(NA, n, n)
diag(dist.mat) <- 0
for(i in 1:(n-1)){
for(j in (i+1):n){
nz <- sum((abs(profiles[, i]) + abs(profiles[, j]))==0)
dc <- sum(abs(profiles[, i] - profiles[, j]) /
(abs(profiles[, i]) + abs(profiles[, j])), na.rm = T) * (k/(k-nz))
dist.mat[i, j] <- dist.mat[j, i] <- dc
}
}
colnames(dist.mat) <- rownames(dist.mat) <- colnames(mat)
dist.mat %>%
knitr::kable(caption = "Table 1.8. Canberra distance calculated using the column profiles of 'seabirds.csv' dataset and the definition given in the exercise 1.E.3: d(p,q) = (k/(k-nz)*(sum(abs(p - q) / (abs(p) + abs(q))))")
(dist_canb <- dist(t(profiles), method = "canberra", diag = T, upper = T)) %>% as.matrix() %>%
knitr::kable(caption = "Table 1.9. Canberra distance calculated using the column profiles of 'seabirds.csv' dataset and the 'dist()' function specifying `method='canberra'`")
| CH | PLI | CI | NS | CL | CT | SI | SPI | SGI | |
|---|---|---|---|---|---|---|---|---|---|
| CH | 0.000000 | 12.73702 | 1.763481 | 20.11905 | 15.09013 | 16.46305 | 22.66969 | 20.077278 | 20.907037 |
| PLI | 12.737016 | 0.00000 | 12.479210 | 19.96096 | 17.37143 | 17.18550 | 22.73001 | 19.751149 | 20.628619 |
| CI | 1.763481 | 12.47921 | 0.000000 | 20.18226 | 14.68790 | 16.91378 | 22.73255 | 20.252762 | 20.990719 |
| NS | 20.119050 | 19.96096 | 20.182256 | 0.00000 | 16.21635 | 19.13141 | 22.65913 | 19.302132 | 20.556911 |
| CL | 15.090126 | 17.37143 | 14.687898 | 16.21635 | 0.00000 | 14.33450 | 22.40510 | 18.434965 | 19.431285 |
| CT | 16.463047 | 17.18550 | 16.913776 | 19.13141 | 14.33450 | 0.00000 | 21.19503 | 19.193085 | 19.299703 |
| SI | 22.669691 | 22.73001 | 22.732545 | 22.65913 | 22.40510 | 21.19503 | 0.00000 | 22.077655 | 21.748062 |
| SPI | 20.077278 | 19.75115 | 20.252762 | 19.30213 | 18.43497 | 19.19308 | 22.07765 | 0.000000 | 9.782611 |
| SGI | 20.907037 | 20.62862 | 20.990719 | 20.55691 | 19.43129 | 19.29970 | 21.74806 | 9.782611 | 0.000000 |
| CH | PLI | CI | NS | CL | CT | SI | SPI | SGI | |
|---|---|---|---|---|---|---|---|---|---|
| CH | 0.000000 | 12.73702 | 1.763481 | 20.11905 | 15.09013 | 16.46305 | 22.66969 | 20.077278 | 20.907037 |
| PLI | 12.737016 | 0.00000 | 12.479210 | 19.96096 | 17.37143 | 17.18550 | 22.73001 | 19.751149 | 20.628619 |
| CI | 1.763481 | 12.47921 | 0.000000 | 20.18226 | 14.68790 | 16.91378 | 22.73255 | 20.252762 | 20.990719 |
| NS | 20.119050 | 19.96096 | 20.182256 | 0.00000 | 16.21635 | 19.13141 | 22.65913 | 19.302132 | 20.556911 |
| CL | 15.090126 | 17.37143 | 14.687898 | 16.21635 | 0.00000 | 14.33450 | 22.40510 | 18.434965 | 19.431285 |
| CT | 16.463047 | 17.18550 | 16.913776 | 19.13141 | 14.33450 | 0.00000 | 21.19503 | 19.193085 | 19.299703 |
| SI | 22.669691 | 22.73001 | 22.732545 | 22.65913 | 22.40510 | 21.19503 | 0.00000 | 22.077655 | 21.748062 |
| SPI | 20.077278 | 19.75115 | 20.252762 | 19.30213 | 18.43497 | 19.19308 | 22.07765 | 0.000000 | 9.782611 |
| SGI | 20.907037 | 20.62862 | 20.990719 | 20.55691 | 19.43129 | 19.29970 | 21.74806 | 9.782611 | 0.000000 |
Realizar un MDS con la distancia de Canberra de R. Comprobar que se trata de una distancia euclídea. Dibujar el mapa. Comparar el resultado con el obtenido con la distancia ji-cuadrado. Utilizar la función procrustes() del paquete vegan.
mds2 <- cmdscale(d = dist_canb, k = 8, eig = T)
mds2$points %>% knitr::kable(caption = "Table 1.10. Coordinates of the MDS from Canberra distances between column profiles shown in table 1.9.")
mds2$eig %>% knitr::kable(caption = "Table 1.11. Eigenvalues of the MDS from Canberra distances between column profiles shown in the table 1.9.")
| CH | -9.004001 | -0.4119965 | 3.403355 | -1.3822228 | -2.9394413 | 2.7220433 | 0.2941079 | 0.6508961 |
| PLI | -6.548694 | -0.5054878 | 3.377353 | -3.1878164 | 7.5606813 | -3.7432170 | -0.4103127 | 0.0303465 |
| CI | -9.199724 | -0.4306546 | 3.359156 | -1.6881892 | -3.2046800 | 1.7311915 | -0.1130372 | -0.6748997 |
| NS | 2.405435 | -0.5506840 | -11.880669 | -6.1205145 | 1.2539361 | 2.3449580 | -0.3070384 | -0.0131285 |
| CL | -2.708016 | -1.5589356 | -5.787390 | 3.8178859 | -4.6417387 | -5.7230153 | -0.5299328 | 0.0591010 |
| CT | -1.427421 | 1.0272123 | -2.678522 | 9.8718390 | 3.9359064 | 2.9963500 | 0.5743877 | -0.0481728 |
| SI | 7.209240 | 14.9136546 | 2.600906 | -0.9028339 | -1.0232724 | -0.6245509 | 0.2122803 | 0.0041246 |
| SPI | 9.087292 | -6.5746806 | 3.163752 | -0.8795986 | -0.5112311 | -0.4970329 | 4.8545370 | -0.0253470 |
| SGI | 10.185890 | -5.9084277 | 4.442059 | 0.4714505 | -0.4301603 | 0.7932733 | -4.5749918 | 0.0170798 |
| x |
|---|
| 462.0537185 |
| 304.9524572 |
| 252.5978579 |
| 166.2239203 |
| 116.1771250 |
| 72.9144505 |
| 45.5147980 |
| 0.8870133 |
| 0.0000000 |
mds2$points %>% as.data.frame() %>% mutate(names = rownames(c$points)) %>%
ggplot(aes(V1, V2)) +
geom_point() +
geom_text(aes(label = names),
position = position_dodge(width = 1),
vjust = -1, hjust = 1, size = 4) +
ggtitle("MDS of Canberra distances",
"of column profiles") +
xlab("1st principal coordinate") + ylab("2nd principal coordinate") +
theme_bw(base_size = 15)
Figure 1.4. MDS of Canberra distances of column profiles. 2D representation (dimensions 1 and 2)
plot(procrustes(c$points, mds2$points, scale = T))
Figure 1.5. MDS of Canberra distances of column profiles. 2D representation (dimensions 1 and 2)
Dos características que tienen las matrices de distancias no Euclídeas es que algunos de los valores propios del escalado multidimensional seran negativos, mientras que algunas de las coordenadas pueden ser numeros complejos. En nuestro caso, ni las coordenadas (tabla 1.10) son numeros complejos, ni los valores propios (tabla 1.11) son negativos. Por este motivo podemos sugerir que las distancias de la tabla 1.9 son Euclídeas.
Realizar un análisis de conglomerados jerárquico con el método deWard2 de la distancia de Canberra según R. Dibujar el dendograma resultante. Dibujar también un heatmap de este análisis con la función heatmap(). Para ello, hay que elegir bien los parámetros distfun= y hclustfun= y una escala de colores. ¿Para qué sirve el heatmap?
hclust_1g <- hcut(x = dist_canb, hc_method = "ward.D2", k = 6)
fviz_dend(hclust_1g, cex = 2 , k = 6, main = "Hierarchical clustering with Canberra distances", sub = "6 clusters",
ggtheme = theme_bw(base_size = 17), show_labels = T, rect = TRUE)
Figure 1.6.
Con el análisis de conglomerados jerárquico (figura 1.6), observamos como muchas de las afirmaciones que hicimos en análisis anteriores salen reforzadas. Para analizar el dendograma de la figura 1.6, debemos observar a que altura se separan cada conglomerado, hasta llegar a los conglomerados (clusters) donde se encuentran las columnas individuales. Cuanto a más altura se de una separción, más lejos estaran los clusteres resultatnes y, para comparar dos clusters, podemos sumar las alturas que los separan:
Más o menos, estudiando el dendograma manualmente (algo que no es eficaç cuando hay muchas columnas), podemos establecer un orden de cercanía entre columnas/conglomerados.
Además, para visualizar un análisis de conglomerados jerárquico también podemos usar un mapa de calor o heatmap. Para esto, podemos usar la función heatmap() del paquete stats, así como heatmap.2() (gplots), pheatmap() (pheatmap) o Heatmap() (ComplexHeatmap)
gplots::heatmap.2(hclust_1g$merge %>% as.matrix())
Figure 1.7.
El siguiente paso es estudiar por algún criterio el número óptimo de conglomerados para el análisis jerárquico. Con la distancia de Canberra según R en particular, lo más sencillo es utilizar el criterio de las siluetas.
Las siluetas son un método para comprobar la consistencia de los conglomerados, de manera que representan como cada observación se parece a su cluster y al cluster más cercano. El valor de la silueta va de -1 a 1, con 1 siendo el valor que indica que más se parece a su propio cluster y -1 el que menos.
Cuando comparamos distintos numeros de clusters, podemos usar el valor medio de la silueta de cada numero de conglomerados para así elegir el número óptimo.
avsi <- c(0)
for(i in 2:8){
hclust_1g <- hcut(x = dist_canb, k = i, isdiss = T, hc_method = "ward.D2")
avsi[i-1] <- hclust_1g$silinfo$avg.width
}
plot(2:8,avsi,type="o", xlab="Number of clusters", ylab="Average silhouette length", main = "Hierarchical clustering")
Figure 1.8.
Como se observa en la figura 1.8, hay dos picos en 4 clústeres (pico local) y en 6 clústers (pico absoluto). Con esto podemos decir que el número òptimo de conglomerados es 6.
Estudiar con la misma distancia el número óptimo de conglomerados con el método PAM.
avsi <- c(0)
for(i in 2:8){
km_vars <- pam(dist_canb, k=i, diss = T)
si <- silhouette(km_vars, dist_canb)
avsi[i-1] <- mean(si[,"sil_width"])
}
plot(2:8,avsi,type="o",xlab="Number of clusters", ylab="Average silhouette length", main = "PAM (K-medioids)")
Figure 1.9.
Curiosamente, usando el método PAM( (figura 1.9), podmeos ver como el pico de la amplitud media de la silueta en los 6 conglomerados desaparace, convirtiéndose en un “valle”, mientras que el pico en los 4 clústeres sigue estando y el pico absoluto está en los 7 grupos.
En este caso, el número óptimo sería los de 7 grupos, aunqué en ambos, los 4 conglomerados parece un número razonable.
avsi <- c(0)
for(i in 2:8){
km_vars <- pam(dist_canb, k=i, diss = T)
si <- silhouette(km_vars, dist_canb)
avsi[i-1] <- mean(si[,"sil_width"])
}
plot(2:8,avsi,type="o",xlab="Number of clusters", ylab="Average silhouette length", main = "PAM (K-medioids)")
Figure 1.10.
De los apartados anteriores se deduce que hay un número razonable de conglomerados, aunque no sea óptimo. Dibujar el dendograma del apartado (g) con esa partición.
hclust_1j <- hcut(x = dist_canb, k = 4, hc_method = "ward.D2", graph = T)
fviz_dend(hclust_1j, cex = 2 , k = 4, main = "Hierarchical clustering with Canberra distances", sub = "4 clusters",
ggtheme = theme_bw(base_size = 17), show_labels = T, rect = TRUE)
Figure 1.11.
Un estudio contiene dos medidas de los anillos de crecimiento en la escala del salmón de Alaska y de Canadá. Los datos se pueden obtener en el libro de Johnson et al. y se adjuntan en el archivo salmon.txt.
salmon <- read.delim("salmon.txt", sep = " ") %>% mutate(Gender = factor(Gender))
salmon_canada <- salmon %>% filter(Origin == "Canadian")
salmon_alaska <- salmon %>% filter(Origin == "Alaskan")
salmon_list <- list(salmon, salmon_canada, salmon_alaska) %>% set_names(nm = c("All salmons", "Canadian salmons", "Alaskan salmons"))
salmon_list_num <- salmon_list %>% purrr::map(~dplyr::select(.data = .x, Freshwater, Marine))
Realizar una estadística descriptiva univariante y multivariante según el factor Origin. Añadir algunos gráficos ilustrativos, en particular el de dispersión.
summary() nos ayuda a obtener información sobre como están estructurados los valores de los datos (centralidad, etc). En este caso vamos a hacerlo en la totalidad de los datos, así como separados por origen.purrr::map(salmon_list, summary)
## $`All salmons`
## Gender Freshwater Marine Origin
## 1:52 Min. : 53.0 Min. :301.0 Alaskan :50
## 2:48 1st Qu.: 99.0 1st Qu.:367.0 Canadian:50
## Median :117.5 Median :396.5
## Mean :117.9 Mean :398.1
## 3rd Qu.:140.0 3rd Qu.:428.2
## Max. :179.0 Max. :511.0
##
## $`Canadian salmons`
## Gender Freshwater Marine Origin
## 1:26 Min. : 90.0 Min. :301.0 Alaskan : 0
## 2:24 1st Qu.:124.0 1st Qu.:345.2 Canadian:50
## Median :140.0 Median :369.5
## Mean :137.5 Mean :366.6
## 3rd Qu.:151.5 3rd Qu.:385.8
## Max. :179.0 Max. :438.0
##
## $`Alaskan salmons`
## Gender Freshwater Marine Origin
## 1:26 Min. : 53.00 Min. :355.0 Alaskan :50
## 2:24 1st Qu.: 86.25 1st Qu.:402.0 Canadian: 0
## Median : 99.00 Median :427.5
## Mean : 98.38 Mean :429.7
## 3rd Qu.:109.00 3rd Qu.:452.5
## Max. :131.00 Max. :511.0
salmon %>% group_by(Origin) %>% summarise(`Mean | Freshwater` = mean(Freshwater),
`Mean | Marine` = mean(Marine),
`SD | Freshwater` = sd(Freshwater),
`SD | Marine` = sd(Marine),
`Median | Freshwater` = median(Freshwater),
`Median | Marine` = median(Marine))
## # A tibble: 2 x 7
## Origin `Mean | Freshwa… `Mean | Marine` `SD | Freshwate… `SD | Marine`
## <fct> <dbl> <dbl> <dbl> <dbl>
## 1 Alask… 98.4 430. 16.1 37.4
## 2 Canad… 137. 367. 18.1 29.9
## # … with 2 more variables: `Median | Freshwater` <dbl>, `Median | Marine` <dbl>
Algo que es fácil de ver es que los salmones Alaskan suelen ser mayores que los Canadian cuando son marinos, pero menores en tamaño cuando se trata de salmones de agua dulce.
var_alaska <- salmon %>% filter(Origin == "Alaskan") %>% dplyr::select(-Gender, -Origin) %>% var()
var_canada <- salmon %>% filter(Origin == "Canadian") %>% dplyr::select(-Gender, -Origin) %>% var()
var_alaska %>% knitr::kable(caption = "Variance-Covariance matrix for Alaskan salmons")
var_canada %>% knitr::kable(caption = "Variance-Covariance matrix for Canadian salmons")
| Freshwater | Marine | |
|---|---|---|
| Freshwater | 260.6078 | -188.0927 |
| Marine | -188.0927 | 1399.0861 |
| Freshwater | Marine | |
|---|---|---|
| Freshwater | 326.0902 | 133.5049 |
| Marine | 133.5049 | 893.2608 |
paste("Varianza total de los salmones de Alaska ", sum(eigen(var_alaska)$values) %>% round(2))
paste("Varianza total de los salmones de Canada ", sum(eigen(var_canada)$values) %>% round(2))
## [1] "Varianza total de los salmones de Alaska 1659.69"
## [1] "Varianza total de los salmones de Canada 1219.35"
paste("Varianza generalizada de los salmones de Alaska ", det(var_alaska) %>% round(2))
paste("Varianza generalizada de los salmones de Canada ", det(var_canada) %>% round(2))
## [1] "Varianza generalizada de los salmones de Alaska 329233.85"
## [1] "Varianza generalizada de los salmones de Canada 273460.04"
paste("Coeficiente de dependencia global de los salmones de Alaska ",
1-det(cor(salmon %>% filter(Origin == "Alaskan") %>% dplyr::select(-Gender, -Origin) )) %>% round(2))
paste("Coeficiente de dependencia global de los salmones de Canada ",
1-det(cor(salmon %>% filter(Origin == "Canadian") %>% dplyr::select(-Gender, -Origin) )) %>% round(2))
## [1] "Coeficiente de dependencia global de los salmones de Alaska 0.1"
## [1] "Coeficiente de dependencia global de los salmones de Canada 0.0600000000000001"
car::scatterplot(salmon_alaska$Freshwater, salmon_alaska$Marine,
main = "Alaskan Freshwater Salmons vs Alaskan Marine Salmons", xlab="Freshwater", ylab = "Marine")
Figure 2.1. Alaskan Freshwater Salmons vs Alaskan Marine Salmons
car::scatterplot(salmon_canada$Freshwater, salmon_canada$Marine,
main = "Canadian Freshwater Salmons vs Canadian Marine Salmons", xlab="Freshwater", ylab = "Marine")
Figure 2.2. Canadian Freshwater Salmons vs Canadian Marine Salmons
car::scatterplot(salmon_alaska$Freshwater, salmon_canada$Freshwater,
main = "Alaskan Freshwater Salmons vs Canadian Freshwater Salmons", xlab="Alaskan", ylab = "Canadian")
Figure 2.3. Alaskan Freshwater Salmons vs Canadian Freshwater Salmons
car::scatterplot(salmon_alaska$Marine, salmon_canada$Marine,
main = "Alaskan Freshwater Marine vs Canadian Marine Salmons", xlab="Alaskan", ylab = "Canadian")
Figure 2.4. Alaskan Freshwater Marine vs Canadian Marine Salmons
Vemos como ni los salmones del mismo lugar (ex. Alaska o Canada) tienen correlación entre agua dulce y agua marina (figuras 2.1 y 2.2), así como los del mismo tipo de agua tienen correlación entre los distintos lugares de origen (figuras 2.3 y 2.4).
salmon %>% melt() %>%
ggplot(aes(Origin, value, color = Origin)) + geom_boxplot() + facet_wrap(~variable, scales = "free") +
theme_pubr(legend = "none") +
ylab("Growth") + xlab("Origin")
Figure 2.4. Boxplot of growth of freshwater and Marine salmons from Alaska and Canada.
Realizar un análisis discriminante lineal. La función lda() del paquete MASS puede servir.
Un análisis discriminante es un método usado para predecir la probabilidad de pertenecer a un grupo determinado basándonos en una o varias variables predictoras. En parte, es similar a la regresión logística (vista en otro tema) en el sentido que ambos métodos son usados para classificar datos en base a variables discretas (ej. grupo 1 y grupo 2), aunqué la regresión logística se usa en casos de variables discretas de dos clases, mientras que el análisis discriminante puede servir para 2 o más clases. 4 También podríamos usar más de una variable como variable predictora.
Cabe destacar que se necesita un grupo de datos conocidos para “entrenar” el modelo y así poder predecir la clasificación de datos nuevos.
Primero hay que usar la función lda() para calcular el modelo, luego, mediante la función predict(), para predecir las probabilidades de cada observación en pertenecer en un grupo u otro (tabla 2.1), así como los coeficientes del model lineal y la predicción de la clase, normalmente según el criterio de probabilidad > 0.5.
En este caso, la variable que defina los grupos será Origin, aunqué podríamos hacerla con otra variable categorica (ej. Gender) o crear una nueva variable que mezcle ambas.
salmon2 <- salmon %>% mutate(group = paste(Origin, Gender, sep = "_"))
lda_2b <- lda(salmon[,2:3], grouping = salmon$Origin)
lda_2b_pred <- lda_2b %>% predict()
lda_2b_pred$posterior %>% round(2) %>% head %>% knitr::kable(caption = "Table 2.1. Probabilities of belonging to Alaskan or Canadian categories for each of the salmons of the dataset. Note that only the first 6 observations are taken.")
| Alaskan | Canadian |
|---|---|
| 0.43 | 0.57 |
| 0.02 | 0.98 |
| 1.00 | 0.00 |
| 1.00 | 0.00 |
| 0.93 | 0.07 |
| 0.99 | 0.01 |
A cada predicción se le asigna una puntuación o score relativo al modelo lineal que se relaciona con la probabilidad de pertenecer a uno u otro grupo y que nos puede dar información sobre la calidad de las predicciones.
ggplot(salmon, aes(lda_2b_pred$x, fill = Origin)) + geom_density(alpha = 0.5) + xlab("scores")
Figure 2.12. Density of the prediction scores in each group (Alaska / Canada)
Como se observa en la figura 2.12, hay un punto del gráfico de densidades en las que los valores de los dos grupos se solapan (cerca del 0), lo que nos puede sugerir que esas predicciones pueden no ser fiables, aunque en su mayor parte, los scores parece buenos.
Otra manera de mirar la calidad del modelo es usando una tabla de clasifiación cruzada, mostrando cuantos de cuandos son predecidos como un grupo son acertados y cuantos no.
table(lda=lda_2b_pred$class,real=salmon$Origin) %>%
knitr::kable(caption = "Table 2.2. Table of cross-classification showing the predictied classes and the real ones.")
| Alaskan | Canadian | |
|---|---|---|
| Alaskan | 44 | 1 |
| Canadian | 6 | 49 |
Vemos en la tabla 2.2 que, de los 50 salmones canadienses, 49 han sido predichos como tal, mientras que en el caso de los salmones de Alaska han sido predichos correctamente en 44 casos, mostrando un poco menos de tasa de acierto, pero igualmente con un resultado aceptable.
Clasificar una observación con un Freshwater de 120 y un valor de Marine de 380.
lda_2b <- lda(formula = Origin ~ Freshwater + Marine, data = salmon)
predict(lda_2b,newdata=data.frame(Freshwater=120,Marine=380))
## $class
## [1] Canadian
## Levels: Alaskan Canadian
##
## $posterior
## Alaskan Canadian
## 1 0.2298261 0.7701739
##
## $x
## LD1
## 1 0.4199577
La predicción es que el salmón es canadiense, mientras que los salmones de medidas parecidas en los datos originales son tanto de Alaska como de Canadá.
Comparar las matrices de covarianzas de las dos poblaciones con el test de la razón de verosimilitudes. También se puede aplicar el test M de Box. Ambos son muy sensibles a la no normalidad de los datos y tienden a rechazar la igualdad de covarianzas.
El análisis discriminante tiene asunciones, una de las cuales es la igualdad entre las mastrices de covarianzas entre las poblaciones que se están comparando. Para comprobar esto, podemos usar un test de verosimilitudes o un test M de Box.
Como sabemos que ambos test son sensibles a la no normalidad de los datos, primero haremos un test de shapiro multivariante para comprobar si estos son normales
mshapiro_test(data = salmon[,2:3])
## # A tibble: 1 x 2
## statistic p.value
## <dbl> <dbl>
## 1 0.995 0.970
Vemos como el p-valor del test de Shapiro-Wilk multivariante es muy proximo a uno, por lo que no podemos rechazar la hipótesis nula de una distribución normal, por lo que podemos fiarnos del resultado de los test.
Test de verosimilitudes:
Test M de Box: el test M de Box, o prueba de box permite comparar las matrices de covarianzas entre dos grupos para ver si estas son iguales. Es un test muy sensible y debido a esto, se recomienda un límite de significancia bastante bajo (i.e. 0.001). También hay que tener en cuenta, que en caso de falta de normalidad de los datos, el resultado del test de Box puede ser significativo debido a esta falta de normalidad y no a la diferencia entre las matrices de covarianzas [explicación sacada de la respuesta de mi PEC1]. En nuestro caso, la distribución parece ser normal, por lo que no debe preocuparnos esta parte. Para llevar a cabo el test, se puede usar la función boxM() del paquete heplots.
boxM_res <- heplots::boxM(salmon[,2:3], salmon$Origin, )
summary(boxM_res)
## Summary for Box's M-test of Equality of Covariance Matrices
##
## Chi-Sq: 10.69615
## df: 3
## p-value: 0.01349
##
## log of Covariance determinants:
## Alaskan Canadian pooled
## 12.70452 12.51891 12.72333
##
## Eigenvalues:
## Alaskan Canadian pooled
## 1 1429.3568 923.1148 1147.0461
## 2 230.3371 296.2362 292.4764
##
## Statistics based on eigenvalues:
## Alaskan Canadian pooled
## product 329233.8474 273460.0441 335483.8619
## sum 1659.6939 1219.3510 1439.5224
## precision 198.3702 224.2669 233.0522
## max 1429.3568 923.1148 1147.0461
Como se puede observar, tanto con la prueba de razón de verosimilitudes que con la prueba de Box, obtenemos que las matrices de covarianzas de los dos grupos son significativamente distintas bajo un nivel de significancia del 95% (p valor inferior a 0.05) y como la muestra parece seguir una distribución normal, esta diferencia es fiable.
Sin embargo, el test M de Box es muy sensible, por lo que se recomienda un límite de significancia de 0.001, por lo que, como el p-valor no es inferior a este, no podemos rechazar con total seguridad la hipótesis nula de covarianzas iguales y, por tanto, el resultado del análisis discriminante lineal es fiable.
En el caso de poblaciones normales con diferentes matrices de covarianzas se clasificará cada observación en el grupo con máxima probabilidad a posteriori, pero entonces las funciones discriminantes no son lineales, ya que tienen un término de segundo grado. Realizar un análisis discriminante cuadrático. La función qda() del paquete MASS nos ayudará
qda_2e <- qda(formula = Origin ~ Freshwater + Marine, data = salmon)
qda_2e_pred <- predict(lda_2b)
table(qda=qda_2e_pred$class,real=salmon$Origin) %>%
knitr::kable(caption = "Table 2.3. Table of cross-classification showing the predictied classes using the quadratic discriminant analysis and the real ones.")
| Alaskan | Canadian | |
|---|---|---|
| Alaskan | 44 | 1 |
| Canadian | 6 | 49 |
En la tabla 2.3 se observan las predicciones de las observaciones originales y se comparan con el grupo original de dichas observaciones. El resultado es el mismo que en la tabla 2.2 pese a que se trata de un análisis discriminante lineal, mientras que en este caso estamos observando el resultado de un análsisi discriminante cuadrático.
Calcular el número de parámetros que hay que estimar en la discriminación lineal y en la cuadrática.
El número de parámetros efectivos a estimar en el LDA se puede calcular con la siguiente función:
\[ K_p + (K-1) \]
Donde \(K\) es el numero de clases de la variable predictora.
Calcular los errores de clasificación con ambas reglas utilizando validación cruzada. Si son similares, nos quedaremos con el análisis lineal que además es más robusto y de mejor interpretación.
Si vemos la tabla 2.2, podemos ver la tabla de datos predichos originales del análisis lineal, mientras que la tabla 2.3 muestra la misma para el análisis cuadrático.
table_lda <- table(lda=lda_2b_pred$class,real=salmon$Origin)
table_lda %>% knitr::kable(caption = "Table 2.2 (rep)")
paste0("Percentage of coincidences without cross-validation: ",round(100*sum(diag(table_lda))/sum(table_lda),2),"%")
paste0("Error rate without cross-validation:: ",round(100*(1-sum(diag(table_lda))/sum(table_lda)),2),"%")
| Alaskan | Canadian | |
|---|---|---|
| Alaskan | 44 | 1 |
| Canadian | 6 | 49 |
## [1] "Percentage of coincidences without cross-validation: 93%"
## [1] "Error rate without cross-validation:: 7%"
table_qda <- table(qda=qda_2e_pred$class,real=salmon$Origin)
table_qda %>% knitr::kable(caption = "Table 2.3. (rep)")
paste0("Percentage of coincidences without cross-validation: ",round(100*sum(diag(table_qda))/sum(table_qda),2),"%")
paste0("Error rate without cross-validation:: ",round(100*(1-sum(diag(table_qda))/sum(table_qda)),2),"%")
| Alaskan | Canadian | |
|---|---|---|
| Alaskan | 44 | 1 |
| Canadian | 6 | 49 |
## [1] "Percentage of coincidences without cross-validation: 93%"
## [1] "Error rate without cross-validation:: 7%"
Vemos como ambas tablas son iguales.
La validación cruzada es un método que nos permite evaluar el resultado de un análisis estadístico y que consiste en repetir el análisis mediante la evaluación de diferentes particiones (grupos de los datos originales que nos permiten entrenar el modelo).
Para llevar a cabo la validación cruzada podemos hacerlo mediante las funciones lda() y qda() indicando el argumento CV = T.
lda_2b_CV <- lda(salmon[,2:3], grouping = salmon$Origin, CV = T)
table_lda_cv <- table(prediction = lda_2b_CV$class, real=salmon$Origin)
table_lda_cv %>% knitr::kable(caption = "Table 2.4. Table of predicted vs original data using LDA with CV.")
paste0("Percentage of coincidences without cross-validation: ",round(100*sum(diag(table_lda_cv))/sum(table_lda_cv),2),"%")
paste0("Error rate without cross-validation:: ",round(100*(1-sum(diag(table_lda_cv))/sum(table_lda_cv)),2),"%")
| Alaskan | Canadian | |
|---|---|---|
| Alaskan | 44 | 1 |
| Canadian | 6 | 49 |
## [1] "Percentage of coincidences without cross-validation: 93%"
## [1] "Error rate without cross-validation:: 7%"
qda_2h_CV <- qda(salmon[,2:3], grouping = salmon$Origin, CV = T)
table_qda_cv <- table(prediction = qda_2h_CV$class, real=salmon$Origin)
table_qda_cv %>% knitr::kable(caption = "Table 2.5. Table of predicted vs original data using QDA with CV.")
paste0("Percentage of coincidences without cross-validation: ",round(100*sum(diag(table_qda_cv))/sum(table_qda_cv),2),"%")
paste0("Error rate without cross-validation:: ",round(100*(1-sum(diag(table_qda_cv))/sum(table_qda_cv)),2),"%")
| Alaskan | Canadian | |
|---|---|---|
| Alaskan | 45 | 3 |
| Canadian | 5 | 47 |
## [1] "Percentage of coincidences without cross-validation: 92%"
## [1] "Error rate without cross-validation:: 8%"
Vemos como, después de usar la validación curzada (CV), los valores de acierto y error de los datos predecidos son muy similares a antes de usar la validación cruzada. Sin embargo, antes de la CV ambos métodos dieron el mismo resultado en cuanto a tasa de acierto y error en las predicciones, pero después pese a ser muy similares, el método LDA parece tener una mejor actuación, por lo que nos quedamos con este.
Alboukadel Kassambara (2017). Correspondence Analysis: Theory and Practice. STHDA. Online aquí.
Michael Greenacre (2008). La pràctica del anàlisis de correspondencias. Capítulo 4: Distancia ji-cuadrado e inercia. Fundación BBVA. Online aquí.
Alboukadel Kassambara (2017). CA in R Using FactoMineR: Quick Scripts and Videos. STHDA. Online aquí.
Alboukadel Kassambara (2018). Discriminant Analysis Essentials in R. STHDA. Online aquí.